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Abstract 

Solvation in l-ethyl-3-methylmidazolium chloride and in l-etliyl-3-methylimidazolium hexaflu- 
orophosphate near equilibrium is investigated via molecular dynamics computer simulations with 
diatomic and benzenelike molecules employed as probe solutes. It is found that electrostriction 
plays an important role in both solvation structure and free energetics. The angular and radial 
distributions of cations and anions become more structured and their densities near the solute 
become enhanced as the solute charge separation grows. Due to the enhancement in structural 
rigidity induced by electrostriction, the force constant associated with solvent configuration fluc- 
tuations relevant to charge shift and transfer processes is also found to increase. The effective 
polarity and reorganization free energies of these ionic liquids are analyzed and compared with 
those of highly polar acetonitrile. Their screening behavior of electric charges is also investigated. 
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I. INTRODUCTION 



Molecular ionic liquids based on bulky cations such as N, N'-dialkylimidazolium and N- 
alkylpyridinium have been the subject of intensive experimental study.— This class of systems 
can exist as a liquid at or near room temperature and is resistant to water. Their vapor 
pressure is essentially zero, so that volatile products can be separated completely through 
distillation and the ionic liquids can be recycled. Accordingly, they provide an environmen- 
tally benign ("green") alternative to toxic organic solvents, which has potentially a wide 
range of applications in organic synthesis and separation chemistry.— 

Structure and solvation properties of these room-temperature ionic liquids (RTILs) have 
been studied via a variety of spectroscopic methods.— I'^'i^i^i^^i^'^i'^^i'^^i'^^i^^ Many solvatochromic 
measurements show that the effective polarity of RTILs is comparable to that of highly 
dipolar solvents such as acetonitrile and small alcohols. Another interesting finding is that 
dynamic Stokes shifts in imidazolium-based RTILs determined with the picosecond time 
resolution are considerably smaller than the corresponding static shifts . Though some- 
what controversial,— this suggests that a substantial part of solvent relaxation in RTILs 
occur on a subpicosecond time scale despite their high viscosity.— 

There have been growing efforts to study RTILs via molecular dynamics (MD) simula- 
tion methods . '^'''i^^i^^i^'^i^^i^^i^'^i^^i^^ Earlier works have focused mainly on pure liquids, such 
as structure and transport properties of RTILs . Recently, we examined equilib- 
rium solvation properties of l-ethyl-3-methylimidazolium chloride (EMI"'"C1~) and 1-ethyl- 
3-methylimidazolium hexafluorophosphate (EMI"'"PFg ) in the presence of a model diatomic 
solute.— We found that in accord with experiments,—!^ the effective polarity of these RTILs 
is very high and nearly comparable to that of ambient water. Solvent fluctuation dynamics 
are characterized by at least two totally different time scales: rapid subpicosecond dynam- 
ics followed by slow non-exponential relaxation due to ion transport.— Furthermore, the 
fast subpicosecond dynamics account for more than 50 % of the entire relaxation of solvent 
fluctuations. Within the context of linear response, this finding lends strong support to 
the existence of ultrafast nonequilibrium solvent relaxation in RTILs based on the indirect 
experimental evidence.i2*i^ 

In this article and its companion^i - hereafter, referred to as II - we extend our previous 
MD study^"^ to investigate solvation structure and dynamics in RTILs both in and out of 
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equilibrium. By employing diatomic and benzenelike solutes as probe molecules, we make 
a detailed analysis of the solvent structural variations with the solute charge distributions 
and their influence on solvation free energetics. We also examine free energy curves along 
a solvent coordinate relevant to charge shift and transfer processes in solution and quantify 
outer-sphere reorganization free energy. Time-dependent Stokes shifts subsequent to an 
instantaneous change in the solute charge distribution are analyzed and comparison is made 
with equilibrium solvent fluctuation dynamics. 

The outline of this paper is as follows: In Sec. II we give a brief account of the model 
description and simulation methods employed in our study. The MD results for solvation 
structure and free energy fluctuations in EMI+Cl" and in EMI+PFg are presented in Sec. Ill, 
while Sec. IV concludes. Dynamic solvation properties are studied in II. 

II. MODEL DESCRIPTION AND SIMULATION METHODS 

The MD simulations are performed using the DL_POLY program.— The simulation cell 
is comprised of a single rigid solute molecule immersed in either EMI+Cl^ or EMI"'"PF^ 
solvent, consisting of 112 pairs of rigid cations and anions. Two different types of solutes 
are considered: a diatomic molecule and a benzenelike one (see Table H}. Two constituent 
atoms of the first type are fixed at a separation of 3.5 A and interact with the solvent through 
Lennard- Jones (LJ) and Coulombic interaction potentials. The LJ parameters, cr = 4A and 
e/^B = 100 K with Boltzmann's constant /cb, and mass, m = 100 amu, are chosen to be 
identical for each constituent atom of the diatomic solute molecule.— Three different charge 
distributions are considered, i.e., a neutral pair (NP) with no charges, an ion pair (IP) with 
unit charge separation (g = ±e), and their intermediate state with g = ±0.5e referred to as 
HIP, where e is the elementary charge. As for the benzenelike solute, we employ cr = 3.5A 
and e/fcB = 40.3K for carbon and o" = 2.5 A and e//cB = 25.2K for hydrogen.— In addition to 
the axially symmetric charge distribution appropriate for the regular ground-state benzene 
(RB), we also study another charge distribution, referred to as "dipolar" benzene (DB); this 
differs from RB by ±e in charge assignments for two C atoms in the para positions while the 
remaining four carbon and six hydrogen atoms have the same partial charges as in the RB 
solute. This yields the dipole moment difference of 13.5 D between the RB and DB charge 
distributions. 
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The solvent parameters employed are the same as in our prior MD study. To be specific, 
we used the AMBER force field^ for the LJ parameters and the partial charge assignments 
of Ref . 3 

for EMI"*". The united atom representation was employed for the CH2 and CH3 
moieties of the ethyl group of the cation as well as for the methyl group (Fig. Q). Experi- 
mental geometry determined by X-ray diffraction^ was used. As for Cl^, we used a = 4.4 A 
and e/A;B = 50.4K. PFg was described as a united atom with a = 5.6 A and e/k^ = 200K. 
The numerical values for the solute and solvent parameters employed in the present study 
are compiled in Table HI 

The simulations are conducted in the canonical ensemble using the Nose-Hoover extended 
system method.— We study two different densities at temperature T = 400 K for each RTIL 
solvent: p = 1.1 and 1.2g/cm3 for EMI+Cr and p = 1.31 and 1.375 g/cm^ for EMI+PFg". 
Periodic, cubic boundary conditions are employed. Long-range electrostatic interactions 
are computed via the Ewald method with a conducting surrounding medium, resulting in 
essentially no truncation of these interactions. The trajectories are integrated via the Verlet 
leapfrog algorithm combined with the quaternion method for rotations with a time step 
of 2fs. Equilibrium simulations are carried out with 2 ns equilibration, followed by a 4 ns 
trajectory from which averages are computed. 



III. RESULTS AND DISCUSSIONS 

The MD results for EMI+Cl" and EMI+PFg are summarized in TableHIland Figs.EHIIl 
We begin our analysis by considering the equilibrium solvent structure around the solutes. 



A. Equilibrium Solvation Structure 

In Fig. 121 radial distribution functions g{r) of the solvent cations and anions around the 
diatomic solutes in EMI+C1~ are displayed. For convenience, the midpoint of two nitrogen 
atoms of the imidazolium ring is employed as the cation center in the present study. The 
most salient feature in Fig. El is that the solvent structure varies markedly with the solute 
charge distribution. As the magnitude of the charges at both sites of the solute is increased, 
their respective first solvation shells consisting mainly of the solvent ions of the opposite 
charge become more pronounced and localized. For example, the Cl~ distribution around the 
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positive site of IP in Fig. vanishes almost completely between its first and second peaks, 
while it shows no substantive structure in the presence of NP. This reveals the development 
of a well-defined shell of Cl^ around the IP solute. As the extent of the solute charge 
separation increases, the position of the first peak in g{r) moves in and its height grows. 
This results in the solvent density increase close to the solute: For instance, the number 
of Cl~ in the first solvation shell of the solute (+) site, i.e., those located within 5 A from 
the site, increases from 1.8 for NP to 3.1 and 4.1 for HIP and IP, respectively (see also 
Fig. El below). We refer to this structure- making in RTILs with increasing solute charge 
separation and accompanying local solvent density enhancement as electrostriction. Similar 
electrostrictive behaviors of dipolar solvents play an important role in various solvation and 
reaction processes involving charge shift and transfer . ^^i^^i^^i'^-^i^^ The basic characteristics 
of the cation distributions in Figs. Efb) and |2tc) are essentially the same as those of the 
anion distributions in Fig. |21^a). We notice that the cation distributions are somewhat less 
structured than the chloride distribution. This is probably due to the extended structure 
and charge distribution of EMI'^ compared with Cl~. 

For additional insight into electrostriction, we consider the integrated ion number N{r) 

N{r) = Attu / g{r') r'^dr' , (1) 
Jo 

which measures the number of anions (or cations) located within a distance r from the solute. 
Here n is the number density of anions (or cations) and g{r) is their radial distribution around 
the solute atom sites. The results for Cl~ and EMI"*" center in the presence of the NP and 
IP solutes are presented in Fig. El Compared with the NP case, N{r) for Cl~ around the 
(+) site of IP is enhanced considerably up to r ^ 6A, while that around the (— ) site is 
somewhat reduced. The qualitative behavior of EMI"*" center in Fig. Efb) is similar to that of 
anions. The strong electrostatic interactions between the IP solute charges and solvent ions 
are responsible for the N{r) trend in the vicinity of the solute. We, however, notice that this 
trend becomes reversed if we consider a larger region around the solute. To be specific, the 
total number of Cl~ present inside a sphere of radius r ^ 6-9 A around the (+) site of IP is 
smaller than that inside the same volume but in the presence of the NP solute. By contrast, 
the total number of Cl~ inside a similar sphere centered at the (— ) site of IP exceeds that 
in the presence of NP. We attribute the former (latter) to the presence of excess negative 
(positive) chages close to the solute positive (negative) site, which repel (attract) Cl~ ions 
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in the outer region (cf. Fig. QHl below). To see this better, we consider a sphere of radius, 
say 4 A, centered at the (+) site of IP (cf. Fig. |2fa)). If we calculate the total charge inside 
this spherical volume including the solute charges, we obtain — 3.2e.— " The positive solute 
charge at the center attracts Cl^ strongly, so that a small neighborhood surrounding the (+) 
site becomes highly negatively-charged. As a consequence, the immdediate vicinity of the 
negatively-charged region is not easily accessible by Cl~ ions from the outer region and thus 
becomes anion-deficient and cation-rich. This is clearly demonstrated in Fig. |2fa), where the 
chloride radial distribution around the (+) site of IP nearly vanishes between r ^ AA and 
6 A. If we consider a larger volume around the solute (+) site to include the above-mentioned 
cation-rich and anion-deficient region, we find that the total charge enclosed inside changes 
the sign and becomes positive; for example, with a radius of 7 A, the total charge becomes 
+3.2e.^^^ This will then attract Cl~ and repel EMI"^ to create an anion-rich and cation- 
deficient region right outside. This shows that the solute charges exert a strong influence 
on solvent structure by inducing anion-rich and cation-deficient, and anion-deficient and 
cation-rich regions that alternate in the solvent medium. The oscillatory character of the 
solvent charge distributions around the solute is studied in detail in Fig. HU] below. 

The cation and anion distributions with respect to the center of the diatomic solutes in 
EMI"''C1~ are exhibited in Fig. |3 It is noteworthy that g{r) in Fig. HJ^a) is considerably less 
structured than that in Fig. Efa). This indicates that the first solvation shell of anions is 
formed mainly around the positive site of the solute (see below). Though much lesser in 
extent, the cations show a similar trend. While the structural variations of the first solvation 
shells are most dramatic, the modulations of the solvent radial distributions by the solute 
charges are not limited to the short range. g{r) in Figs. El and |3] and N{r) in Fig. El vary 
substantially with the solute electronic structure for r ~ 10 A. This implies that the influence 
of the solute charges on the solvent ions extends well beyond a separation of ~ 10 A. We will 
come back to this point below. 

Figure shows the angular probability distribution P(cos6') of EMI"*" and of Cl^ around 
the diatomic solute, where 9 denotes the angle between the solute dipole vector and the 
solute-to-ion center-to-center direction. (For example, 6 = 0° and 180° correspond to the 
directions of the positive and negative sites of the solute from its center, respectively.) The 
ions are grouped together in three different regions according to the distance r from the 
solute center: region I (r < S.SA), region II (5.5 < r < qA) and region III (r > qA). 
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In each region, the probabihty is normahzed to unity. In the presence of the NP solute, 
P{cos6) for both cations and anions in regions II and III are found to be independent of 6; 
i.e., their distributions are isotropic. In region I, on the other hand, while nearly constant 
around cos^ = 0, P{cos6) vanishes as cos^^ approaches ±1. This anisotropy arises from the 
presence of the solute diatoms at 6* = 0° and 180°. Because of their extended structure, 
the cation distribution in region I around cos 6' = is not as isotropic as the corresponding 
anion distribution. 

The angular distributions of the ions around the IP solute in Figs. El^c) and El^d) differ 
markedly from those in Figs. Efa) andEfb). The former are much more structured than 
the latter as in the case of the radial distribution functions. The narrow width of P{cos6) 
in region I in Figs. E^c) and Efd) confirms that the anions and cations close to the solute 
(cf. Fig. are situated mainly around the positive and negative sites of the IP solute, 
respectively. These ions constitute primarily the first solvation shells corresponding to the 
first main peaks of the anion and cation radial distributions in Fig. |2l In region II, P{cos6) 
for both cations and anions are rich in structure with multiple peaks, manifesting that their 
distributions are highly anisotropic. It is interesting to note that the cation and anion 
populations in regions I and II tend to fluctuate in cos 6* in a staggered fashion, viz., they 
oscillate out of phase. This indicates that EMI"*" rings and Cl~ ions tend to alternate in their 
angular distributions around the IP solute for r < 9 A. The ionic distributions in region III 
are considerably more isotropic than those in regions I and II. However, the former are not 
as isotropic as their counterpart in the presence of NP. This is another indication that the 
effects of the solute charge on the solvent structure goes considerably beyond a separation 
of ~ 10 A (see below). 

The MD results for the solvent structure around the diatomic solutes in EMI+PFg are 
shown in Fig.lHl As expected, EMI+PFg shows essentially the same electrostrictive behavior 
as EMI+Cl". We nonetheless notice that the first peaks of the anion distributions around the 
IP and HIP are located at larger values of r in the former RTIL than in the latter because 
PFg is bigger than Cl^. In addition, the anion peak heights are reduced in EMI+PFg 
compared with EMI+Cl^. By contrast, the cation distributions show the opposite trend. To 
be specific, g{r) for the Ml site and cation center around the IP and HIP solutes become 
enhanced in EMI+PFg at both p = 1.31 and 1.375 g/cm'^, compared with those in EMI+Cl" 
at p = 1.1 and 1.2g/cm'^. Though not presented here, the distribution of the C2 site (Fig.[TI) 
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shows a similar trend. This imphes that in the density range we investigated, the cation 
contributions to RTIL energetics and dynamics would be more significant in EMI+PFg than 
in EMI+Cl". We will revisit this point in Sec. 1111 Bl below and in II. We also note that the 
behavior of N{r), Eq. in EMI+PFg" (not shown here) is similar to that in EMI+Cl^ in 
Fig. El 

We now turn to benzenelike solutes. The cation and anion radial distributions in the 
presence of RB and DB in EMI+Cl" and in EMI+PFjT are exhibited in Fig. [7| Comparison 
with Figs. 121 and (HI shows that their basic trends are very similar to the diatomic solute 
case. Specifically, both the anion and cation distributions become more structured and their 
densities very close to the solute become higher as the solute charge separation is increased. 
Also as the anion size grows, the cation structure near the solute becomes enhanced, es- 
pecially in the case of DB. One noteworthy feature is that in both RTILs, the ratio of the 
first peak heights of the EMI"*" center and anion radial distributions is larger in the pres- 
ence of the benzenelike DB than that in the presence of the diatomic IP. For example, in 
EMI+C1~, the main anion peak around IP in Fig. Efa) is about 4-5 times higher than the 
corresponding cation peak in Fig. |21^c), while the former [Fig. jHl^a)] is only twice as large as 
the latter [Fig. Etb)] in the presence of DB. In the case of EMI+PFg , the first peak of EMI"*" 
is actually higher than that of PFg in the presence of DB, whereas it is the opposite in the 
presence of IP. Thus, the relative enhancement of cation structure with increasing solute 
charge separation is more pronounced near the benzenelike solute than near the diatomic 
solute. This implies that cations play a more significant role in solvation in the presence of 
the former solute than in the presence of the latter (see below). 

We examine the ring orientation of EMI"*" relative to benzene as a function of the angle 
between the two unit vectors normal to the benzene and imidazolium rings. The results for 
EMI+PFg are displayed in Fig. |H1 While the ring orientation of EMI"*" in the "bulk" region, 
i.e., more than 5.5 A from the solute center, is almost isotropic, those in the first solvation 
shell show a rather anisotropic distribution with peaks around cos0 = ±0.9 for RB and 
±0.98 for DB. This means that relatively parallel configurations of the EMI"*" and benzene 
rings, i.e., 7r-stacking, are favored near the solute, consonant with a previous simulation 
study. Though not presented here, our analysis shows that on the average, the Ml and El 
groups of EMI"*" (Fig. ^ close to the solute face outward when viewed from the center of 
the solute. This allows a closer approach of the cations to the solute by avoiding the steric 
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hindrance between the two alkyl groups and benzene ring. 

Finally, we investigate the screening behavior of electric charges in EMI+C1~ and 
EMI"'"PFg . In the absence of a rigorous theoretical framework to describe RTILs, we em- 
ploy the formulation developed for simple ionic systems. In the strongly-coupled regime, the 
charge distribution around ionic species a at a large separation takes the forn>SLSLS 

ga(r) = ^e-^/^in(^ + ^) , (2) 

where Aa describes the a-dependence, A is the screening length, and d and ip are the period 
and phase shift associated with charge oscillations in the r-space, respectively. We examine 
Qa{r) with respect to both the anionic species and various charge sites of EMI+. For 
illustration, the MD results for ln|r(5Q!(r)| for a = Cl^ and C2 site of the cation (Fig. 
in EMI+C1~ at density p = l.lg/cm'^ are plotted in Fig. |H1 Despite the limitation in the 
r-range probed via MD in the present study, we observe two notable features there. First, to 
a good approximation, the r-dependence of the envelop of \rQa{r)\ displays an exponential 
decay except at short separations. In fact, for r > 6A, Eq. (0) apparently provides quite 
an accurate description for the charge distributions in RTILs, including oscillations (see 
below). The respective screening lengths for a = Cl~ and C2 are around 9 and 10 A. A 
similar analysis for the Ml and Nl sites of EMI"*" yields A ~ 7 and 9 A, respectively. We 
thus conclude that the screening length in EMI+C1~ at p = l.lg/cm'^ is in the range 7- 
10 A. For EMI+PF^ at 1.31 g/cm^, we obtain A = 10-20 A. Though numerical uncertainties 
associated with the latter are somewhat large, our results undoubtedly demonstrate that the 
screening length is of the order of 10 A for both RTILs. We also notice that electric charges 
are better screened in EMI+Cl^ with small anions than in EMI+PFg with large anions. 
For perspective, we make brief contact with the well-known Debye-Hiickel theory,'^^*^ which 
yields 

Qa{r) = ^e-^/'-- (3) 
r 

with the screening length Adh = {'^'^Yla^o'^a)~^^'BT. Its direct application to the RTILs 
would yield Adh ~ 0.15 A and 0.18 A for EMI+Cr and EMI+PFg , respectively. The Debye- 
Hiickel theory, which completely ignores molecular size effects and is thus valid only in the 
infinite dilution limit of ionic solutions, overestimates the ionic screening behavior of RTILs 
by about two orders of magnitude.— 

The second aspect we point out is that as alluded to above, the charge distribution 
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oscillates quite regularly in Fig. El The respective d values for EMI+Cl" and EMI"'"PF^ 
are about 5.7 A and 6.6 A. This seems to indicate that the distribution of concentric (but 
smeared-out) shells of alternating charges around a central ion, a picture originally put 
forward for molten salts, also applies to the present RTILs. While this may be the case 
for spherical anions, we expect that charge distributions around nonspherical EMI"^ with an 
extended molecular structure would be rather aniosotropic. 

The radial charge distributions Q±{r) around the positive and negative sites of the IP 
solute in EMI+Cl^ are exhibited in Fig. ^| Analogous to the charge distributions around 
the solvent ions in Fig.lHl Q±{r) show oscillations in r with a period c? ~ 6 A. This alternation 
of cation-rich and anion-rich regions around IP was already noted above in connection with 
N{r) in Fig. |21 We notice that the locations of the minima coincide with the peak 

positions of the chloride distribution around the (+) site of IP in Fig. Efa). The locations 
of the Q-{r) maxima generally agree with the peak positions of the cation distributions 
in Figs. Efb) and|2fc).— For comparison, the solvent charge distribution in the presence 
of NP is also displayed in Fig. E3 In direct contrast with the IP case, the radial charge 
density around NP is essentially neutral with minor fluctuations. This difference in Q{r) 
between the NP and IP emphasizes the dramatic effect the solute charges have on solvent 
structure and related charge distributions in ionic liquids. We also notice that the IP charges 
are screened by the solvent just like the ions. The screening lengths for Q±{r) appear to 
be similar to those associated with the solvent ions, viz., ~ 10 A (although we would need 
considerably improved MD statistics to accurately quantify the former). This provides an 
explanation for the relatively long range of the solute charge effect, i.e., the structure of 
the solvent ions located more than 10 A from the solute is still significantly affected by the 
solute charge distributions. We note that the screening of the solute charges is invoked to 
explain the non-fiuctuating behavior of the average solute-solvent electrostatic interactions 

n 

at large separations in Ref. |22l 
B. Solvation Free Energetics 

We proceed to free energetic aspects of solvation in RTILs, especially equilibrium solvation 
stabilization and fluctuations relevant to charge transfer and shift processes. For simplicity, 
we assume that the solute is characterized by two nonpolarizable electronic states a and fe. 
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which are degenerate in energy. For a given solvent configuration Q and a solute electronic 
state i (= a, b), we denote the total energy of the combined solute-solvent system as Ei{Q). 
The Franck-Condon (FC) energy AEa^b{Q) associated with the a ^ b transition in the 
presence of Q is given by 

AEa-.b{Q) = Ea{Q) - Eb{Q) . (4) 

Since Eq. ^ measures the solvation stabilization difference between states a and b, which 
varies with the solvent configuration Q, the FC energy AEa-^b provides a very convenient 
variable ("solvent coordinate") that can be used in lieu of Q to describe the collective 
infiuence of the solvent on the solute.— Since all intermolecular interactions are assumed to 
be pairwise additive and the solute LJ parameters do not vary with the solute electronic 
states, only the Coulombic interactions between the solute and solvent contribute to AEa^b 
in Eq. ^ in our description. 

The probability distribution Pa/b of the FC energy AEa^b when the solvent is in equilib- 
rium with the a-state solute is given by 

Pa/biAE,^b) = jdQ f:\Q) K^Ea^b - AEa^biQ)) . (5) 

Here /^"^(Q) is the equilibrium ensemble distribution of Q in the presence of the a-state 
solute and S{- ■ ■) is the Dirac delta function. The equilibrium average and fiuctuations of 
AEa^biQ) then read 

(AEa^b) = J dxxPa/b{x) 
{{5AEa^,Y) = j dxx'Pa/bix) - {AEa^,)^ , (6) 

where 5AEa^b is the deviation of AEa^b from its equilibrium average, i.e., 5AEa^b = 
AEa-*b — {AEa^b) ■ We also introduce a free-energy function for solvent fiuctuations: 

Fa/biAEa^,) = -ksT In P,/fe( A^,^,) . (7) 

This defines an effective electronic potential energy curve, upon which solvation processes 
of the a-state solute, monitored via a FC transition to the 6-state, occur. With the aid of 
the equipartition principle, we can determine the force constant ka/b associated with Fa/b 
according to^ 

^^1' - WAE-bf) ■ 
11 



The underlying assumption in Eq. (jH)) is that Pa/bi^Ea^b) defined in Eq. @ is gaussian and 
thus Fa/b{^Ea-,b) is harmonic. 

Since (AEa^b) is the (free) energy difference between solute state a in equilibrium with 
the solvent and its FC-(de)excited state b and the two states are degenerate in energy in 
vacuum, it describes the energy shift in the a ^ b transition induced by solvation in RTILs. 
Therefore, (AEa^b) with dipolar a-state is similar to various empirical polarity scales based 
on solvent spectral shifts, which gauge the solvation power of RTILs. To quantify the effective 
polarity of the RTILs, we have computed (AEip^^p) and (Ai?DB^RB) for EMI+C1~ and for 
EMI^PFg and compiled the results, together with those for acetonitrile for comparison, in 
Table ITTl It is found that the effective polarity of these RTILs is very high. In fact, their 
values for (AEip^-^p) are larger than that of highly dipolar acetonitrile, indicating that as 
far as solvation stabilization is concerned, EMI+C1~ and EMI'^'PFg' behave as more "polar" 
solvents than CH3CN. This is in good agreement with recent spectroscopic measurements^"^ 
because l-butyl-3-methylimidazolium hexafiuorophosphate (BMI+PFg ) and EMI"'"PFg are 
similar in their effective polarity (see below). We also observe that the values of (Ai^DB^RB) 
are smaller than the corresponding (AEip^-^p) values by more than 10 %. We believe that 
the difference in the degree of charge separation between the diatomic and benzenelike 
solutes is mainly responsible for this. Since the dipole moment difference of the DB and 
RB solutes (13.5 D) is smaller than that of the IP and NP solutes (16. 7D), the associated 
solvation stabilization difference is smaller for the former solute type than the latter. We 
also notice that the value of (Ai^ip^Np) and thus effective polarity increase with the solvent 
density. 

Another noteworthy feature about (AEa^b) is that its values for IP NP and DB RB 
are lower in EMI+PFg than in EMI"''C1~. This suggests that for a given cationic species, 
the solvating power of RTILs decreases as their anion size grows. To gain further insight, 
we have analyzed the cation and anion contributions, denoted by (AE^^f^) and (AE^^^i^), 
to (AEa-^b)- The results in Table HTl show that the anions play a major role in solvation 
stabilization in EMI+Cl^. For both a = IP and DB, {AE^^ arising from Cl~ is larger 
than {AE'^f^ from EMI+ by more than a factor of two. In EMI+PFg", on the other hand, 
the anion contributions become lowered significantly by 30-40 % for a = IP and by ~ 50 % 
for DB, compared with EMI"*"C1~. This has the origin in the fact that bulky PFg cannot 
approach the (+) site of the solute as closely as small Cl~ [cf. Figs. Efa) andlH^a) above]. By 
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contrast, the cation contributions to (AEa^b) increase with the anion size. We attribute this 
to the above-mentioned enhancement in cation structure near the solute, especially around 
its negative site, when Cl~ is replaced by PFg [see Figs. I21b)(c) and|ni^b)(c) above]. As 
a consequence, the contributions of cations to solvation stabilization and thus to effective 
solvent polarity become comparable to (and sometimes even exceed) those of anions in 
EMI"'"PFg . Nevertheless, effects of common cations on (AEa^b) changes are of secondary 
importance compared with those of the anions which we vary. By analogy we would expect 
that modulations of the cation size by, e.g., varying its alkyl chain length would have similar 
effects on the solvation strength and effective polarity of RTILa^ though they would be 
considerably smaller than those we found here between EMI+Cl" and EMI"''PFg . Thus, for 
example, the effective polarity of BMI"^PFg with a longer butyl chain would be somewhat 
lower than that of EMI+PFg with a shorter ethyl chain. 

We notice in Table |n] that for a given RTIL, the relative contribution of cations to 
(AEa^b) compared to that of anions is larger for DB RB (i.e., benzenelike solute) than 
for IP NP (i.e., diatomic solute). In EMI+Cl", for example, {AE^%) is ~29 % of the total 
spectral shift for IP NP, while it accounts for 32 % in the case of DB RB. In EMI+PFg" 
at 1.375 g/cm'^, the cation contributions for the former and latter are, respectively, 46% and 
56 %. We believe that this is closely related to differing solvent structure enhancement with 
the solute type as discussed near Fig. El above. To be specific, for the RTILs studied here, 
the relative degree of structure making associated with cations is higher near the benzenelike 
solute than near the diatomic solute as the solute charge separation grows. Accordingly we 
expect that the relative contribution to (AEa^b) from cations is larger in the presence of 
DB than that in the presence of IP. This suggests that respective roles played by cations and 
anions in (AEa^b) and related solvation properties are influenced by the solute type — both 
geometry and charge distribution — mainly through the modulations of the ion distributions 
near the solute. This would mean that the relative contributions of anions and cations to the 
solvent spectral shift in the presence of, e.g., betaine-30 studied in R.ef. I27L are reduced and 
enhanced, respectively, compared to the diatomic solute considered here. This is because 
the main positive site of the former is more embedded in the solute and thus less exposed 
to the solvent ions than its negative site. 

Turning to solvent fluctuations relevant to solvation in RTILs, we have found that 
{{SAEa^bY) diminishes as the solute charge separation is increased f Table IIT|). This means 
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that the effective potential energy curve Fa/b{^Ea^b) in Eq. ((Tj) becomes tighter and its 
force constant ka/b in Eq. (jH)) larger with the growing solute charge separation (see Fig. [TT|l. 
This trend is more pronounced in the presence of benzenelike solutes than diatomic solutes. 
We ascribe this to the enhancement of rigidity in solvent structure close to the solute. As 
mentioned several times above, electrostriction yields solvent density enhancement near the 
solute as the charge separation of the latter increases. We believe that this reduces the free 
volume available for individual solvent ions and thus suppresses solvent structural fluctua- 
tions near the solute. This in turn lowers {{5AEa^bY) ■, which is directly related to solvent 
configuration fluctuations. It is interesting to note that according to the present and pre- 
vious studies,—!^ aprotic solvents such as CH3CN and CH3CI show a similar trend in the 
solvent force constant. 

We add a cautionary remark: The solvent density increase in the proximity of the solute 
also strengthens the solute-solvent Coulombic interactions, i.e., the magnitude of /S.Ea^h 
increases. Assuming all other aspects of solvation would remain unchanged, we can infer 
that this enhancement in the electrostatic interaction would increase {{SAEa^bY) ■ This is 
the opposite of the structural rigidity effects mentioned above. Thus qualitatively speaking, 
electrostriction causes two opposing factors to increase with the solute charge separation — 
viz., solvent structural rigidity and electrostatic interaction strength, which play antagonistic 
roles in the fluctuations of AEa^b- For RTILs (and acetonitrile) studied here, the modulation 
of structural rigidity appears to dominate over that of electrostatic interaction strength, so 
that {{SAEa^bY) decreases as the solute charge separation grows. We parenthetically note 
that the opposite trend obtains in water, i.e., {{SAEa-^bY) increases with the solute charge 
separation . ^^1'^'^ 

We next consider the solvent reorganization free energy \a/b, which describes the free- 
energy cost associated with rearranging the solvent on the Fa/b curve from the configurations 
in equilibrium with solute state a to those with solute state b. Since the solvent force constant 
ka/b varies with the solute electronic state, so does Xa/b- If we assume that Fa/b is harmonic, 
the state-dependent reorganization free energy \a/b is given by^ 

Xa/, = \ka/, (A^n + Ai^r-ia)' ^ \ka/, {{AEa^b) + {AEb^a)f , (9) 

where Ai?^^ is the value of AEa^b at the minimum of Fa/b^ Using the local force constant 
determined via Eq. (jH}, we calculate Xa/b in Eq. (jH)) and present the results in Table ITTl For 
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EMI+Cl" at p = l.lg/cm^, we obtain Ajp/np = 66.8 kcal/mol and Anp/ip = 48.5 kcal/mol. 
The corresponding values for EMI"^PF^ at 1.375 g/cm^ are 68.4 and 45. 7 kcal/mol. This 
shows that charge shift processes in RTILs such as IP ^ NP are accompanied by huge sol- 
vent reorganization. Here it is of interest to make comparison with acetonitrile, in which 
the reorganization free energy is Ajp/np = 41.5 and Anp/ip = 36.9 kcal/mol at room temper- 
ature. According to recent spectroscopic measurements,—!^ the reorganization free energy in 
BMI+PFg is comparable to that in acetonitrile. Since BMI"'"PFg is expected to be somewhat 
lower in effective polarity than EMI+PFg (see above), the solvent reorganization energy in 
the former would be smaller than that in the latter [see Eq. (jl(J|) below]. This is consistent 
with our finding that \a/b in EMI+Cl" and EMI"'"PFg are in general larger than that in 
CH3CN. We also mention that because the free energy curve is tighter in the presence of IP 
than in the presence of NP, Aip/np is larger than Anp/ip and similarly for the DB and RB 
states.— 

In view of the state-dependent solvation behavior observed above, we briefly examine the 
assumption of harmonic potentials invoked in Eq. (jU)). The reorganization free energies and 
FC transition energies satisfj*^ 

K/b + Aft/a = + A^ria ~ {^Ea^b) + (AEft^„), (10) 

whether or not Fa/b and Fb/a are harmonic. Comparing the predictions for \a/b + \/a in 
the harmonic assumption with the MD results for (AEa^b) + {^Eb^a) we find easily that 
Eq. (fTUj) holds only approximately. For instance, the estimate for (AE'ip^np) + (Ai^NP^ip) 
in the harmonic approximation is 115kcal/mol in EMI+Cl^ at p = l.lg/cm'^, while the 
MD yields lOlkcal/mol. The corresponding results for EMI'*"PFg at 1.375g/cm^ are 114 
and 95 kcal/mol, respectively. We thus conclude that the anharmonicity in Fa/b plays a 
non-negligible role in solvation free energetics. In view of the drastic solvent structure 
change and related solvent force constant variation with the solute charge distribution, 
the presence of nonlinearity in solvent response is not surprising. What is surprising is 
that despite huge electrostrictive effects, the harmonic approximation for Fa/b with the 
state-dependent force constant^ provides a reasonable description to quantitate the solvent 
reorganization free energy in RTILs. While a similar harmonic approximation works well 
for dipolar solvents,—!^ the RTILs are a totally different matter in that they are comprised 
of free ions. Related issues for solvent djTiamics, e.g., linear response, are considered in II. 
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Before we conclude, we consider the united-atom description employed in the present 
study for perspective. Because the CH3 and CH2 moieties of EMI"^ (Fig- are spherically 
symmetric in their shape and charge distribution in the united-atom description, the direc- 
tionality in their interactions with other molecules (i.e., dependence on relative orientations 
of their H atoms) present in real ionic liquids is completely neglected in our MD. Also two 
rigid alkyl groups of EMI"*" could introduce spurious hindrance compared with a more ac- 
curate flexible description. As for PFg , while the rigid, united-atom representaion may be 
reasonable for short-range interactions because of its high symmetry, the extended charac- 
ter of its charge distribution ignored in the present study could play a nonnegligible role in, 
e.g., /S.Ea^b- To gain insight into thses issues, we have performed a test simulation using 



the flexible, all-atom description of Ref. 
Ref. 



22I for EMI"*" and the rigid, all-atom description of 



23|for PFg . Its results, averaged over a 2 ns trajectory in the presence of the diatomic 
solute in EMI"*" PFg are given in Table ITTl We notice that the united-atom and all-atom de- 
scriptions compare well in their predictions for solvation properties. In fact, considering the 
difference in their parametrization, the degree of agreement between the two in the solvent 
spectral shift and its anion and cation components, solvent fluctuations and reorganization 
free energy is quite remarkable. Though very preliminary, we have found a similar agreement 
between the united-atom and all-atom descriptions for EMI"^C1~. This seems to indicate 
that despite its approximate nature, the united-atom description employed in the present 
study captures the essential features of RTILs very well. For related issues on dynamics, the 
reader is referred to II. 



IV. CONCLUDING REMARKS 



We have studied equilibrium solvation structure and free energetics in ionic liquids con- 
sisting of l-ethyl-3-methylimidazolium cations paired with either chloride or hexafluorophos- 
phate anions. We have found that solvent structure varies dramatically with the solute 
charge distribution. To be specific, as the solute charge separation increases, both the cation 
and anion radial distribution functions become more structured and their densities close to 
the solute become enhanced. Also their angular distributions around the solute become 
highly anisotropic except in the bulk region. For instance, the cation and anion populations 
near the IP solute, though smeared-out, tend to alternate in angular distributions. The 
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charge distributions around the solvent ions and IP solute show screened oscillations. The 
screening length is estimated to be of the order of 10 A in both EMI~^C1~ and EMI"*'PFg 
although the statistical error associated with the screening of the IP solute is large. It was 
observed that the screening length is not sensitive to the presence of a solute or its charge 
character. 

We examined the effective polarity of EMI+Cl^ and EMI+PFj^, measured as solvent 
spectral shifts that gauge the solvation stabilization of dipolar solutes. The MD results show 
that the effective polarity of EMI+C1~ and EMI+PFfT is very high; insofar as their solvating 
strengths are concerned, both RTILs behave as a more polar solvent than acetonitrile, in 
accord with experiments.—*^ EMI+Cl" is more polar of the two, indicating that solvation 
power decreases with the increasing ion size. Related solvent reorganization free energies 
were found to be in general higher in these RTILs than in acetonitrile. 

We have analyzed respective roles played by the anions and cations in solvation. While 
both of them make important contributions to solvent polarity and reorganization free en- 
ergy, their relative importance tends to vary with the solute type. To be specific, the relative 
contribution of anions becomes reduced in the presence of the benzenelike solute, compared 
with the diatomic solute case. In parallel with this, the relative enhancement in cation 
structure near the solute tends to be more pronounced in the presence of the former than 
in the presence of the latter solute. This suggests that variations in solute structure and 
charge distributions have a considerable influence on details of solvation by modulating the 
ion distributions close to the solute. 

To gain insight into the united-atom description employed in the present study, we con- 
ducted a test simulation with an all-atom description for EMI"'"PFg and compared with the 
united-atom results. We obtained good agreement between the two in the solvent spectral 
shift, its anion and cation contributions, equilibrium solvent fluctuations and reorganization 
free energy (related issues on solvent dynamics are examined in II). Very preliminary results 
for EMI+Cl" in the all- atom representation show a similar agreement. This indicates that 
the united-atom description employed here provides a decent framework to simulate RTILs. 

The free-energy curves relevant to solvent fluctuations and associated force constants 
were found to vary with the solute charge distribution. To be specific, the free-energy curves 
become tighter as the solute charge separation grows. This was attributed to the rigidity 
enhancement in solvent structure induced by electrostriction. This demonstrates that a 
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linear solvation scheme is in general not valid for RTILs. It was found that the harmonic 
approximation with state-dependent force constants works reasonably well to quantify, e.g., 
the solvent reorganization free energy.— 

Considering the huge electrostrictive effects on the RTIL structure observed in this work, 
we expect that dynamic properties of solute molecules, e.g., rotational and vibrational energy 
relaxation times, will vary significantly with the solute charge distribution. For example, 
the rotational friction of the solute will grow considerably with its dipole moment analogous 
to dielectric friction in dipolar solvents. A similar trend is also expected for the solute 
vibrational energy relaxation. An analysis of these properties will be reported elsewhere.— 

In the following paper,— we consider equilibrium and nonequilibrium solvation dynamics 
in RTILs. 
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TABLE I: Lennard-Jones parameters and partial charges 





atom 


(Tii (A) 


en (kJ/mol) 


Qi (e) 


EMI+ 


N 


3.250 


0.7113 


-0.2670 




C2 


3.400 


0.3598 


0.4070 




N 


3.250 


0.7113 


-0.2670 




C4 


3.400 


0.3598 


0.1050 




C5 


3.400 


0.3598 


0.1050 




Ml 


3.905 


0.7330 


0.3160 




El 


3.800 


0.4943 


0.2400 




M3 


3.800 


0.7540 


0.0760 




H2 


2.421 


0.0628 


0.0970 




H4 


2.421 


0.0628 


0.0940 




H5 


2.421 


0.0628 


0.0940 


cr 




4.400 


0.4190 


-1.000 






5.600 


1.6680 


-1.000 


diatomic 


NP 


4.000 


0.8298 


0.0000 




IP 


4.000 


0.8298 


±1.0000 


benzene*^ 


C 


3.500 


0.3352 


-0.1350 




H 


2.500 


0.2094 


0.1350 



The LJ parameters are from Ref. iio! and partial charges are from Ref. 
Ref. 0. 



2G 
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TABLE II: MD results'') 



solvent 


density 


a/b 




i^Eltb) 


i^Eltb) 


{{5^Ea-.b?) 


A /I 

a/b 


EMI+Cl- 


1.1 


NP/IP 


0.25 


-0.04 


0.29 


83.3 


48.5 






IP/NP 


100.6 


71.1 


29.5 


60.5 


66.8 


EMi+cr 


1.1 


RB/DB 


-0.41 


-0.69 


0.28 


77.5 


38.4 






DB/RB 


86.9 


59.3 


27.6 


54.2 


54.9 


EMi+cr 


1.2 


NP/IP 


-0.5 


-1.65 


1.15 


72.3 


56.1 






IP/NP 


101.5 


71.7 


29.8 


57.9 


70.0 


EMI+PFg 


1.31 


NP/IP 


0.20 


0.43 


-0.23 


79.8 


35.1 






IP/NP 


83.8 


43.7 


40.1 


64.5 


43.5 


EMI+PFg 


1.34 


NP/IP 


-0.44 


-1.02 


0.58 


80.9 


37.3 






IP/NP 


87.6 


48.3 


39.3 


61.2 


49.3 


EMI+PFg 


1.375 


NP/IP 


0.53 


0.60 


-0.07 


77.8 


45.7 






IP/NP 


94.1 


51.0 


43.1 


52.0 


68.4 


EMI+PFg 


1.375 


RB/DB 


0.15 


0.28 


-0.13 


81.2 


25.5 






DB/RB 


72.0 


31.8 


40.2 


38.7 


53.5 


CHsCN'^) 


0.73 


NP/IP 


0.46 






46.4 


38.0 






IP/NP 


76.5 






41.3 


42.7 



°^ Units for energy: kcal/mol. Due to the difference in MD statistics, some results here 
are somewhat different from those in Ref IitI. 

^) Estimation in the harmonic approximation with Eq. (jHI). 

'^•' Preliminary results obtained with an all-atom description. See the text for details. 
'^^ MD results at T = 300 K. The LJ parameters and partial charges employed are from 
Ref.y. 
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FIG. 1: Structure of the l-ethyl-3-methylimidazolium cation. In simulations, the united atom 
representation is employed for Ml, El and M3 moieties. 

FIG. 2: Radial distributions of cations and anions around the diatomic solute in EMI+Cl^ at 
p = l.lg/cm^: (a) Cl~ around the solute (+) site; (b) Ml around the solute (— ) site; (c) EMI+ 
center around the solute (— ) site. Here r is the solute- ion distance in units of A and the midpoint 
of two N sites in the imidazolium ring (cf. Fig. ^ is employed as the EMI"*" center in (c). The 
results for g{r) at higher density p = 1.2g/cm^ are almost the same as those at p = l.lg/cm^ 
shown here. 

FIG. 3: Number of (a) Cl^ and (b) EMI+ ions surrounding the diatomic solute in EMI+Cl^ at 
p = l.lg/cm^. The dotted and dashed lines represent, respectively, N{r) around the positive and 
negative sites of IP, while the solid line denotes that around the neutral atoms of NP. In (b), the 
EMI+ center is used in the calculations of N{r). 

FIG. 4: Radial distributions of cations and anions with respect to the center of the diatomic solute 
in EMI+Cr at p = l.lg/cm^: (a) CI" and (b) EMI+ center. 

FIG. 5: Angular probability distributions of (a) anions and (b) cations around the NP solute and 
(c) anions and (d) cations around the IP solute in EMI"''C1~ at p = l.lg/cm^. Solid, dashed and 
dotted lines denote the probability distributions in regions I, II, and III, respectively, while 9 is 
the angle between the IP dipole vector and the solute-to-ion center-to-center direction. In (a) and 
(b), P{cos9) for region III is nearly the same as that for region II and thus is not shown. 

FIG. 6: Radial distributions of (a) anions around the (+) site, (b) Ml around the (— ) site and (c) 
EMI"*" center around the (— ) site of the diatomic solute in EMI+PFg at p = 1.375 g/cm'^. Radial 
distributions at lower density p = 1.31 g/cm^ are somewhat less structured but their essential 
features including the peak positions remain nearly the same. 
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FIG. 7: Radial distributions of solvent ions around the carbon sites of the benzenelike solute. The 
results for anions and the center of cations in EMI+Cl^ are presented in (a) and (b), respectively, 
and the corresponding distributions in EMI+PFg are shown in (c) and (d). In the case of DB, the 
anion and cation distributions around only those C atoms that are, respectively, more positive and 
negative in their charge assignments than the remaining four are shown here. The solvent densities 
are the same as in Figs. [U and 1^1 

FIG. 8: Probability distributions of the EMI+ ring orientation with respect to the benzenelike 
solute in EMI+PFg at T = 400 K and p = 1.375 g/cm'^. cj) is the angle between the normal vectors 
of the imidazolium and benzene rings. The solid and dashed lines represent P(cos (j)) in the first 
solvation shells of DB and RB, respectively, which consist of the cations with the center-to-center 
distance of less than 5.5 A from the solute. The rest are defined as the "bulk" region, P(cos 0) of 
which in the presence of RB is plotted as a dotted line. P(cos</>) is normalized to unity in each of 
the two regions. The distribution of the DB ring orientation in the bulk is nearly the same as that 
of RB and thus is not shown here. 

FIG. 9: Absolute value oirQa{r) (a =01", C2) in EMI+Cr at p = l.lg/cm^ as a function of the 
radial distance r in the presence of NP. The dotted lines have the slopes —0.115 and —0.10 for C2 
and CI", respectively, corresponding to the inverse screening length A"^. Though not presented 
here, comparison with other cases — specifically, pure RTILs and RTIL solutions containing IP — 
shows that the screening length and oscillation period are rather insensitive to the presence of the 
diatomic solute and its charge character. 

FIG. 10: Radial charge distributions Q±{r) around the positive ( — ) and negative ( ) site of 

the IP solute in EMI+Cr at p = 1.1 g/cm^. 

FIG. 11: Free energy curves in the presence of (a) a diatomic and (b) a benzenelike solute in 
EMI^^Cl" at p = l.lg/cm'^. Both the free energy F^/fe s-^icl deviation 6AEa^b of AEa^b from its 
average are given in units of kcal/mol. 
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